home *** CD-ROM | disk | FTP | other *** search
/ AGA Toolkit '97 / The AGA Toolkit '97.iso / miscellaneous / science / maths / calc / source / comfunc.c < prev    next >
Encoding:
C/C++ Source or Header  |  1996-09-07  |  11.3 KB  |  596 lines
  1. /*
  2.  * Copyright (c) 1993 David I. Bell
  3.  * Permission is granted to use, distribute, or modify this source,
  4.  * provided that this copyright notice remains intact.
  5.  *
  6.  * Extended precision complex arithmetic non-primitive routines
  7.  */
  8.  
  9. #include "cmath.h"
  10.  
  11.  
  12. /*
  13.  * Round a complex number to the specified number of decimal places.
  14.  * This simply means to round each of the components of the number.
  15.  * Zero decimal places means round to the nearest complex integer.
  16.  */
  17. COMPLEX *
  18. cround(c, places)
  19.     COMPLEX *c;
  20.     long places;
  21. {
  22.     COMPLEX *res;        /* result */
  23.  
  24.     res = comalloc();
  25.     res->real = qround(c->real, places);
  26.     res->imag = qround(c->imag, places);
  27.     return res;
  28. }
  29.  
  30.  
  31. /*
  32.  * Round a complex number to the specified number of binary decimal places.
  33.  * This simply means to round each of the components of the number.
  34.  * Zero binary places means round to the nearest complex integer.
  35.  */
  36. COMPLEX *
  37. cbround(c, places)
  38.     COMPLEX *c;
  39.     long places;
  40. {
  41.     COMPLEX *res;        /* result */
  42.  
  43.     res = comalloc();
  44.     res->real = qbround(c->real, places);
  45.     res->imag = qbround(c->imag, places);
  46.     return res;
  47. }
  48.  
  49.  
  50. /*
  51.  * Compute the result of raising a complex number to an integer power.
  52.  */
  53. COMPLEX *
  54. cpowi(c, q)
  55.     COMPLEX *c;        /* complex number to be raised */
  56.     NUMBER *q;        /* power to raise it to */
  57. {
  58.     COMPLEX *tmp, *res;    /* temporary values */
  59.     long power;        /* power to raise to */
  60.     unsigned long bit;    /* current bit value */
  61.     int sign;
  62.  
  63.     if (qisfrac(q))
  64.         math_error("Raising number to non-integral power");
  65.     if (zge31b(q->num))
  66.         math_error("Raising number to very large power");
  67.     power = (zistiny(q->num) ? z1tol(q->num) : z2tol(q->num));
  68.     if (ciszero(c) && (power == 0))
  69.         math_error("Raising zero to zeroth power");
  70.     sign = 1;
  71.     if (qisneg(q))
  72.         sign = -1;
  73.     /*
  74.      * Handle some low powers specially
  75.      */
  76.     if (power <= 4) {
  77.         switch ((int) (power * sign)) {
  78.             case 0:
  79.                 return clink(&_cone_);
  80.             case 1:
  81.                 return clink(c);
  82.             case -1:
  83.                 return cinv(c);
  84.             case 2:
  85.                 return csquare(c);
  86.             case -2:
  87.                 tmp = csquare(c);
  88.                 res = cinv(tmp);
  89.                 comfree(tmp);
  90.                 return res;
  91.             case 3:
  92.                 tmp = csquare(c);
  93.                 res = cmul(c, tmp);
  94.                 comfree(tmp);
  95.                 return res;
  96.             case 4:
  97.                 tmp = csquare(c);
  98.                 res = csquare(tmp);
  99.                 comfree(tmp);
  100.                 return res;
  101.         }
  102.     }
  103.     /*
  104.      * Compute the power by squaring and multiplying.
  105.      * This uses the left to right method of power raising.
  106.      */
  107.     bit = TOPFULL;
  108.     while ((bit & power) == 0)
  109.         bit >>= 1L;
  110.     bit >>= 1L;
  111.     res = csquare(c);
  112.     if (bit & power) {
  113.         tmp = cmul(res, c);
  114.         comfree(res);
  115.         res = tmp;
  116.     }
  117.     bit >>= 1L;
  118.     while (bit) {
  119.         tmp = csquare(res);
  120.         comfree(res);
  121.         res = tmp;
  122.         if (bit & power) {
  123.             tmp = cmul(res, c);
  124.             comfree(res);
  125.             res = tmp;
  126.         }
  127.         bit >>= 1L;
  128.     }
  129.     if (sign < 0) {
  130.         tmp = cinv(res);
  131.         comfree(res);
  132.         res = tmp;
  133.     }
  134.     return res;
  135. }
  136.  
  137.  
  138. /*
  139.  * Calculate the square root of a complex number, with each component
  140.  * within the specified error.  If the number is a square, then the error
  141.  * is zero.  For sqrt(a + bi), this calculates:
  142.  *    R = sqrt(a^2 + b^2)
  143.  *    U = sqrt((R + abs(a))/2)
  144.  *    V = b/(2 * U)
  145.  *    then sqrt(a + bi) = U + Vi if a >= 0,
  146.  *    or abs(V) + sgn(b) * U  if a < 0
  147.  */
  148. COMPLEX *
  149. csqrt(c, epsilon)
  150.     COMPLEX *c;
  151.     NUMBER *epsilon;
  152. {
  153.     COMPLEX *r;
  154.     NUMBER *A, *B, *R, *U, *V, *tmp1, *tmp2, *epsilon2;
  155.     long m, n;
  156.  
  157.     if (ciszero(c) || cisone(c))
  158.         return clink(c);
  159.     if (cisreal(c)) {
  160.         r = comalloc();
  161.         if (!qisneg(c->real)) {
  162.             r->real = qsqrt(c->real, epsilon);
  163.             return r;
  164.         }
  165.         tmp1 = qneg(c->real);
  166.         r->imag = qsqrt(tmp1, epsilon);
  167.         qfree(tmp1);
  168.         return r;
  169.     }
  170.  
  171.     A = qlink(c->real);
  172.     B = qlink(c->imag);
  173.     n = zhighbit(B->num) - zhighbit(B->den);
  174.     if (!qiszero(A)) {
  175.         m = zhighbit(A->num) - zhighbit(A->den);
  176.         if (m > n)
  177.             n = m;
  178.     }
  179.     epsilon2 = qscale(epsilon, n/2);
  180.     R = qhypot(A, B, epsilon2);
  181.     qfree(epsilon2);
  182.     if (qisneg(A))
  183.         tmp1 = qsub(R, A);
  184.     else
  185.         tmp1 = qadd(R, A);
  186.     qfree(A);
  187.     tmp2 = qscale(tmp1, -1L);
  188.     qfree(tmp1);
  189.     U = qsqrt(tmp2, epsilon);
  190.     qfree(tmp2);
  191.     qfree(R);
  192.     if (qiszero(U)) {
  193.         qfree(B);
  194.         qfree(U);
  195.         return clink(&_czero_);
  196.     }
  197.     tmp1 = qdiv(B, U);
  198.     V = qscale(tmp1, -1L);
  199.     qfree(tmp1);
  200.     r = comalloc();
  201.     if (qisneg(c->real)) {
  202.         if (qisneg(B)) {    
  203.             tmp1 = qneg(U);
  204.             qfree(U);
  205.             U = tmp1;
  206.             tmp2 = qabs(V);
  207.             qfree(V);
  208.             V = tmp2;
  209.         }
  210.         r->real = V;
  211.         r->imag = U;
  212.     } else {
  213.         r->real = U;
  214.         r->imag = V;
  215.     }
  216.     qfree(B);
  217.     return r;
  218. }
  219.  
  220.  
  221. /*
  222.  * Take the Nth root of a complex number, where N is a positive integer.
  223.  * Each component of the result is within the specified error.
  224.  */
  225. COMPLEX *
  226. croot(c, q, epsilon)
  227.     COMPLEX *c;
  228.     NUMBER *q, *epsilon;
  229. {
  230.     COMPLEX *r;
  231.     NUMBER *a2pb2, *root, *tmp1, *tmp2, *epsilon2;
  232.  
  233.     if (qisneg(q) || qiszero(q) || qisfrac(q))
  234.         math_error("Taking bad root of complex number");
  235.     if (cisone(c) || qisone(q))
  236.         return clink(c);
  237.     if (qistwo(q))
  238.         return csqrt(c, epsilon);
  239.     r = comalloc();
  240.     if (cisreal(c) && !qisneg(c->real)) {
  241.         r->real = qroot(c->real, q, epsilon);
  242.         return r;
  243.     }
  244.     /*
  245.      * Calculate the root using the formula:
  246.      *    croot(a + bi, n) =
  247.      *        cpolar(qroot(a^2 + b^2, 2 * n), qatan2(b, a) / n).
  248.      */
  249.     epsilon2 = qscale(epsilon, -8L);
  250.     tmp1 = qsquare(c->real);
  251.     tmp2 = qsquare(c->imag);
  252.     a2pb2 = qadd(tmp1, tmp2);
  253.     qfree(tmp1);
  254.     qfree(tmp2);
  255.     tmp1 = qscale(q, 1L);
  256.     root = qroot(a2pb2, tmp1, epsilon2);
  257.     qfree(a2pb2);
  258.     qfree(tmp1);
  259.     tmp1 = qatan2(c->imag, c->real, epsilon2);
  260.     qfree(epsilon2);
  261.     tmp2 = qdiv(tmp1, q);
  262.     qfree(tmp1);
  263.     r = cpolar(root, tmp2, epsilon);
  264.     qfree(root);
  265.     qfree(tmp2);
  266.     return r;
  267. }
  268.  
  269.  
  270. /*
  271.  * Calculate the complex exponential function to the desired accuracy.
  272.  * We use the formula:
  273.  *    exp(a + bi) = exp(a) * (cos(b) + i * sin(b)).
  274.  */
  275. COMPLEX *
  276. cexp(c, epsilon)
  277.     COMPLEX *c;
  278.     NUMBER *epsilon;
  279. {
  280.     COMPLEX *r;
  281.     NUMBER *tmp1, *tmp2, *epsilon2;
  282.  
  283.     if (ciszero(c))
  284.         return clink(&_cone_);
  285.     r = comalloc();
  286.     if (cisreal(c)) {
  287.         r->real = qexp(c->real, epsilon);
  288.         return r;
  289.     }
  290.     epsilon2 = qscale(epsilon, -2L);
  291.     r->real = qcos(c->imag, epsilon2);
  292.     r->imag = qlegtoleg(r->real, epsilon2, _sinisneg_);
  293.     if (qiszero(c->real)) {
  294.         qfree(epsilon2);
  295.         return r;
  296.     }
  297.     tmp1 = qexp(c->real, epsilon2);
  298.     qfree(epsilon2);
  299.     tmp2 = qmul(r->real, tmp1);
  300.     qfree(r->real);
  301.     r->real = tmp2;
  302.     tmp2 = qmul(r->imag, tmp1);
  303.     qfree(r->imag);
  304.     qfree(tmp1);
  305.     r->imag = tmp2;
  306.     return r;
  307. }
  308.  
  309.  
  310. /*
  311.  * Calculate the natural logarithm of a complex number within the specified
  312.  * error.  We use the formula:
  313.  *    ln(a + bi) = ln(a^2 + b^2) / 2 + i * atan2(b, a).
  314.  */
  315. COMPLEX *
  316. cln(c, epsilon)
  317.     COMPLEX *c;
  318.     NUMBER *epsilon;
  319. {
  320.     COMPLEX *r;
  321.     NUMBER *a2b2, *tmp1, *tmp2;
  322.  
  323.     if (ciszero(c))
  324.         math_error("Logarithm of zero");
  325.     if (cisone(c))
  326.         return clink(&_czero_);
  327.     r = comalloc();
  328.     if (cisreal(c) && !qisneg(c->real)) {
  329.         r->real = qln(c->real, epsilon);
  330.         return r;
  331.     }
  332.     tmp1 = qsquare(c->real);
  333.     tmp2 = qsquare(c->imag);
  334.     a2b2 = qadd(tmp1, tmp2);
  335.     qfree(tmp1);
  336.     qfree(tmp2);
  337.     tmp1 = qln(a2b2, epsilon);
  338.     qfree(a2b2);
  339.     r->real = qscale(tmp1, -1L);
  340.     qfree(tmp1);
  341.     r->imag = qatan2(c->imag, c->real, epsilon);
  342.     return r;
  343. }
  344.  
  345.  
  346. /*
  347.  * Calculate the complex cosine within the specified accuracy.
  348.  * This uses the formula:
  349.  *    cos(a + bi) = cos(a) * cosh(b) - sin(a) * sinh(b) * i.
  350.  */
  351. COMPLEX *
  352. ccos(c, epsilon)
  353.     COMPLEX *c;
  354.     NUMBER *epsilon;
  355. {
  356.     COMPLEX *r;
  357.     NUMBER *cosval, *coshval, *tmp1, *tmp2, *tmp3, *epsilon2;
  358.     int negimag;
  359.  
  360.     if (ciszero(c))
  361.         return clink(&_cone_);
  362.     r = comalloc();
  363.     if (cisreal(c)) {
  364.         r->real = qcos(c->real, epsilon);
  365.         return r;
  366.     }
  367.     if (qiszero(c->real)) {
  368.         r->real = qcosh(c->imag, epsilon);
  369.         return r;
  370.     }
  371.     epsilon2 = qscale(epsilon, -2L);
  372.     coshval = qcosh(c->imag, epsilon2);
  373.     cosval = qcos(c->real, epsilon2);
  374.     negimag = !_sinisneg_;
  375.     if (qisneg(c->imag))
  376.         negimag = !negimag;
  377.     r->real = qmul(cosval, coshval);
  378.     /*
  379.      * Calculate the imaginary part using the formula:
  380.      *    sin(a) * sinh(b) = sqrt((1 - a^2) * (b^2 - 1)).
  381.      */
  382.     tmp1 = qsquare(cosval);
  383.     qfree(cosval);
  384.     tmp2 = qdec(tmp1);
  385.     qfree(tmp1);
  386.     tmp1 = qneg(tmp2);
  387.     qfree(tmp2);
  388.     tmp2 = qsquare(coshval);
  389.     qfree(coshval);
  390.     tmp3 = qdec(tmp2);
  391.     qfree(tmp2);
  392.     tmp2 = qmul(tmp1, tmp3);
  393.     qfree(tmp1);
  394.     qfree(tmp3);
  395.     r->imag = qsqrt(tmp2, epsilon2);
  396.     qfree(tmp2);
  397.     qfree(epsilon2);
  398.     if (negimag) {
  399.         tmp1 = qneg(r->imag);
  400.         qfree(r->imag);
  401.         r->imag = tmp1;
  402.     }
  403.     return r;
  404. }
  405.  
  406.  
  407. /*
  408.  * Calculate the complex sine within the specified accuracy.
  409.  * This uses the formula:
  410.  *    sin(a + bi) = sin(a) * cosh(b) + cos(a) * sinh(b) * i.
  411.  */
  412. COMPLEX *
  413. csin(c, epsilon)
  414.     COMPLEX *c;
  415.     NUMBER *epsilon;
  416. {
  417.     COMPLEX *r;
  418.  
  419.     NUMBER *cosval, *coshval, *tmp1, *tmp2, *epsilon2;
  420.  
  421.     if (ciszero(c))
  422.         return clink(&_czero_);
  423.     r = comalloc();
  424.     if (cisreal(c)) {
  425.         r->real = qsin(c->real, epsilon);
  426.         return r;
  427.     }
  428.     if (qiszero(c->real)) {
  429.         r->imag = qsinh(c->imag, epsilon);
  430.         return r;
  431.     }
  432.     epsilon2 = qscale(epsilon, -2L);
  433.     coshval = qcosh(c->imag, epsilon2);
  434.     cosval = qcos(c->real, epsilon2);
  435.     tmp1 = qlegtoleg(cosval, epsilon2, _sinisneg_);
  436.     r->real = qmul(tmp1, coshval);
  437.     qfree(tmp1);
  438.     tmp1 = qsquare(coshval);
  439.     qfree(coshval);
  440.     tmp2 = qdec(tmp1);
  441.     qfree(tmp1);
  442.     tmp1 = qsqrt(tmp2, epsilon2);
  443.     qfree(tmp2);
  444.     r->imag = qmul(tmp1, cosval);
  445.     qfree(tmp1);
  446.     qfree(cosval);
  447.     if (qisneg(c->imag)) {
  448.         tmp1 = qneg(r->imag);
  449.         qfree(r->imag);
  450.         r->imag = tmp1;
  451.     }
  452.     return r;
  453. }
  454.  
  455.  
  456. /*
  457.  * Convert a number from polar coordinates to normal complex number form
  458.  * within the specified accuracy.  This produces the value:
  459.  *    q1 * cos(q2) + q1 * sin(q2) * i.
  460.  */
  461. COMPLEX *
  462. cpolar(q1, q2, epsilon)
  463.     NUMBER *q1, *q2, *epsilon;
  464. {
  465.     COMPLEX *r;
  466.     NUMBER *tmp, *epsilon2;
  467.     long scale;
  468.  
  469.     r = comalloc();
  470.     if (qiszero(q1) || qiszero(q2)) {
  471.         r->real = qlink(q1);
  472.         return r;
  473.     }
  474.     epsilon2 = epsilon;
  475.     if (!qisunit(q1)) {
  476.         scale = zhighbit(q1->num) - zhighbit(q1->den) + 1;
  477.         if (scale > 0)
  478.             epsilon2 = qscale(epsilon, -scale);
  479.     }
  480.     r->real = qcos(q2, epsilon2);
  481.     r->imag = qlegtoleg(r->real, epsilon2, _sinisneg_);
  482.     if (epsilon2 != epsilon)
  483.         qfree(epsilon2);
  484.     if (qisone(q1))
  485.         return r;
  486.     tmp = qmul(r->real, q1);
  487.     qfree(r->real);
  488.     r->real = tmp;
  489.     tmp = qmul(r->imag, q1);
  490.     qfree(r->imag);
  491.     r->imag = tmp;
  492.     return r;
  493. }
  494.  
  495.  
  496. /*
  497.  * Raise one complex number to the power of another one to within the
  498.  * specified error.
  499.  */
  500. COMPLEX *
  501. cpower(c1, c2, epsilon)
  502.     COMPLEX *c1, *c2;
  503.     NUMBER *epsilon;
  504. {
  505.     COMPLEX *tmp1, *tmp2;
  506.     NUMBER *epsilon2;
  507.  
  508.     if (cisreal(c2) && qisint(c2->real))
  509.         return cpowi(c1, c2->real);
  510.     if (cisone(c1) || ciszero(c1))
  511.         return clink(c1);
  512.     epsilon2 = qscale(epsilon, -4L);
  513.     tmp1 = cln(c1, epsilon2);
  514.     tmp2 = cmul(tmp1, c2);
  515.     comfree(tmp1);
  516.     tmp1 = cexp(tmp2, epsilon);
  517.     comfree(tmp2);
  518.     qfree(epsilon2);
  519.     return tmp1;
  520. }
  521.  
  522.  
  523. /*
  524.  * Return a trivial hash value for a complex number.
  525.  */
  526. HASH
  527. chash(c)
  528.     COMPLEX *c;
  529. {
  530.     HASH hash;
  531.  
  532.     hash = qhash(c->real);
  533.     if (!cisreal(c))
  534.         hash += qhash(c->imag) * 2000029;
  535.     return hash;
  536. }
  537.  
  538.  
  539. /*
  540.  * Print a complex number in the current output mode.
  541.  */
  542. void
  543. comprint(c)
  544.     COMPLEX *c;
  545. {
  546.     NUMBER qtmp;
  547.  
  548.     if (_outmode_ == MODE_FRAC) {
  549.         cprintfr(c);
  550.         return;
  551.     }
  552.     if (!qiszero(c->real) || qiszero(c->imag))
  553.         qprintnum(c->real, MODE_DEFAULT);
  554.     qtmp = c->imag[0];
  555.     if (qiszero(&qtmp))
  556.         return;
  557.     if (!qiszero(c->real) && !qisneg(&qtmp))
  558.         math_chr('+');
  559.     if (qisneg(&qtmp)) {
  560.         math_chr('-');
  561.         qtmp.num.sign = 0;
  562.     }
  563.     qprintnum(&qtmp, MODE_DEFAULT);
  564.     math_chr('i');
  565. }
  566.  
  567.  
  568. /*
  569.  * Print a complex number in rational representation.
  570.  * Example:  2/3-4i/5
  571.  */
  572. void
  573. cprintfr(c)
  574.     COMPLEX *c;
  575. {
  576.     NUMBER *r;
  577.     NUMBER *i;
  578.  
  579.     r = c->real;
  580.     i = c->imag;
  581.     if (!qiszero(r) || qiszero(i))
  582.         qprintfr(r, 0L, FALSE);
  583.     if (qiszero(i))
  584.         return;
  585.     if (!qiszero(r) && !qisneg(i))
  586.         math_chr('+');
  587.     zprintval(i->num, 0L, 0L);
  588.     math_chr('i');
  589.     if (qisfrac(i)) {
  590.         math_chr('/');
  591.         zprintval(i->den, 0L, 0L);
  592.     }
  593. }
  594.  
  595. /* END CODE */
  596.